Force network ensemble: a new approach to static granular matter 
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An ensemble approach for force distributions in static granular packings is developed. This 
framework is based on the separation of packing and force scales, together with an a-priori flat 
measure in the force phase space under the constraints that the contact forces are repulsive and 
balance on every particle. We show how the formalism yields realistic results, both for disordered and 
regular triangular "snooker ball" configurations, and obtain a shear-induced unjamming transition 
of the type proposed recently for athermal media. 
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The fascinating properties of static granular matter 
are closely related to the organization of the interparticle 
contact forces into highly heterogeneous force networks 
0. The probability density of contact forces, P{f), has 
emerged as a key characterization of a wide range of ther- 
mal and athermal systems 10: l^t tSt lE S • Most of 
these studies so far focussed on the broad tail of this dis- 
tribution. Recently, however, the non-universal shoulder 
of P{f) for small forces has received increasing atten- 
tion, since it appears to contain nontrivial information 
on the state of the system: P(/) exhibits a peak at some 
small value of / for "jammed" systems which gives way to 
monotonic behavior above the glass transition 0, H, ■ 
This hints at a possible connection between jamming, 
glassy behavior and force network statistics, and under- 
scores the paramount importance of developing a theoret- 
ical framework for the statistics and spatial organization 
of the forces. 

A first step towards a statistical description of force 
networks is the definition of a suitable ensemble over 
which to take averages. A popular approach is that 
of Edwards, who proposed to take an equal probabil- 
ity for all "blocked" states of a given energy and density 
|lCllll| . Each point in this ensemble defines a contact ge- 
ometry and a configuration of (balancing) contact forces. 
Even though the precise particle locations and contact 
forces are related, the crucial point is that for hard par- 
ticles (most granular matter, hard-sphere colloids, parti- 
cles with a steep Lennard- Jones interactions) a separa- 
tion of scales occurs : forces inside a pile of steel balls 
originate from minute deformations (~ 10^''). 

In this Letter we exploit this scale separation by fo- 
cussing on the ensemble of force configurations for a given 
fixed packing configuration - see Fig. 1. As we will show, 
the forces can be considered underdetermined in this ap- 
proach, since the number of unknown forces exceeds the 
number of balance equations: such packings are refered 
to as hyperstatic. Both packings of entirely rigid, fric- 
tional particles and packings of frictionless particles that 
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FIG. 1: Two different mechanically stable force configurations 
for (a) a "snooker-triangle packing" and for (b) an irregular 
contact network of 1024 particles (only part is shown); the 
thickness of the lines is proportional to the contact force. The 
"force network ensemble" samples all possible force configu- 
rations for a given contact network with an equal probability. 



deform ~ 10^'^ are hyperstatic For transparancy, 

however, we will restrict ourselves to two dimensional 
frictionless hyperstatic packings. In the spirit of Ed- 
wards, we sample all allowed force configurations for a 
given packing with equal probability. Interestingly, the 
idea to restrict the ensemble to fixed geometries has also 
been suggeste d by Bouchaud in the context of extremely 
weak tapping [ISj . 

Our ensemble approach captures the essential features 
of force networks (Fig. 1) and leads to new insights on 
the non-universal "shoulder" of P{f)- In addition, by 
separating the contact geometry from the forces, we can 
start to disentangle the separate roles of contact and 
stress anisotropies in sheared systems. The conceptual 
advantage of not averaging over packing geometries is 
complemented by practical advantages: our protocal is 
numerically cheap and analytically accessible. 

Formulation We study 2D packings of N frictionless 
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disks of radii Ri with centers r^. We denote the interpar- 
ticle force on particle i due to its contact with particle j 
by fij. There are zN/2 contact forces in such packings {z 
being the average contact number), and for purely repul- 
sive central forces we can write fy = fijVij /\rij\, where 
all fij (= fji) are positive scalars. For a fixed contact 
geometry, we are thus left with 2N unknowns and 
zN/2 unknowns fij These satisfy the conditions of 
mechanical equilibrium, 

27V eqs.: ^ = 0, where = r, - r^, (1) 

3 ■' 

and once a force law F is given, the forces are explicit 
functions of the particle locations: 

zN/2 eqs.: /„■ = F{r,f,R„Rj) . (2) 

Pakings of infinitely hard particles have z = 4 and are 
thus isostatic: For rigid particles Eqs. lj2Jl reduce to zN/2 
constraints on the 2N coordinates which can only be 
satisfied if 2 < 4, while Eqs. (Q) can only be solved if 
z > 4; combining these yields z = 4 jlSlllq . 

However, for particles of finite hardness, packings 
are typically hyperstatic with z > 4. A key param- 
eter which quantifies the separation of length scales is 
e = where ( ) denotes an average over the 

packing. We will avoid the strict isostatic e = case, but 
focus instead on the regime where e is small and varia- 
tions of the force of order (/) result in minute variations 
of r,j , of relative size e. Hence, for e <^1, Eqs. (1) and (2) 
can be considered separated, and the essential physics is 
then given by the force balance constraints Eqs. (Q) with 
fixed Tj. In this interpretation, there are more degrees 
of freedom {zN/2) than constraints {2N), leading to an 
ensemble of force networks - see Fig. 1. 

It is important to note that different points in our en- 
semble do not correspond to precisely the same packing 
of exactly the same particles. Our stochastic approach 
describes different force configurations arising in, e.g., 
experiments on "regular" packings of imperfect cannon 
balls or packings under weak tapping Experi- 
mentally, it has become clear that the macroscopic prop- 
erties of granular packings are sensitive to many (coarse 
grained) parameters such as local densities, anisotropics 
and contact numbers, and that it is very difficult to estab- 
lish the relevant characteristics of a packing. Our ensem- 
ble averages only over microscopic variations of the pack- 
ing, which have a strong effect on the local forces but not 
on macroscopic properties, while keeping the important 
characteristics fixed. The restriction to a minute part of 
the Edwards ensemble may therefore help to disentangle 
the separate roles of contact and stress anisotropics. 

The ensemble of force networks for a fixed contact ge- 
ometry is constructed as follows, (i) Assume an a-priori 
fiat measure in the force phase space {/}. (ii) Impose 




FIG. 2: Interparticle force P{f) for various triangular 
"snooker" packings (Fig. la). The inset shows the evolution 
of the tail for large systems. 

the 2A'^ linear constraints given by the mechanical equi- 
librium Eqs. (Hi) Consider repulsive forces only, i.e., 
V/ij > 0. (iv) Set an overall force scale by applying a 
fixed pressure, similar to energy or particle number con- 
straints in the usual thermodynamic ensembles. 

We will first illustrate our formalism for a simple tri- 
angular snooker- like packing (Fig. ^). Then, for an ir- 
regular packing of 1024 particles (Fig. QJd), we compare 
our ensemble approach to MD simulations by varying the 
inverse "hardness" e. The ensemble reproduces the P{f) 
for sufficiently hard particles well. Finally, we find that 
applying a shear stress yields an "unjamming" transition 
in our framework. 

Regular packings The triangular snooker-like pack- 
ings shown in Fig.^ have 3N unknown forces (boundary 
forces included) that are constrained by the 2A^ equations 
of mechanical equilibrium. Even though the packing 
geometry is completely regular, the ensemble approach 
yields irregular force networks and a broad P{f)- 

Labeling each bond by a single index fc, the mechanical 
equilibrium can be expressed as 

A/ = with /= (/i,/2,---,/3Jv) , (3) 

where A is a 3A^ times 2A^ sparse matrix. There is thus 
an A^-dimensional subspace of allowed force configura- 
tions that obey mechanical equilibrium. Eq. is homo- 
geneous, but in physical realizations an overall force scale 
is determined by the externally applied stresses and/or 
the gravitational bulk forces. The simplest manner to 
do so here is to fix the external pressure by specifying 
the total boundary forces. For the snooker-triangles it 
then follows that the sum of all forces is constant. We 
arc thus considering the phase space defined by the force 
balance the "pressure" constraint fk ~ i^tot, and 
the condition that all /'s are positive: 

Af^b and V A>0, (4) 

where the fixed matrix A is the matrix A extended by 
the pressure constraint and b = (0, 0, 0, • ■ • , 0, Ftot). 
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FIG. 3: Comparison of the P{f) obtained by our sampling of 
a frozen geometry, and MD simulations of increasingly hard 
particles under constant pressure in the limit T —> 0. 



To compute P{f) for larger packings, we have applied 
a simulated annealing procedure T^. Starting from an 
ensemble of random initial force configurations we sam- 
ple the space of mechanically stable networks, using a 
penalty function whose degenerate ground states are so- 
lutions of Eq. Q . We have carefully checked that results 
do not depend on the initial configurations, and further- 
more perfectly reproduce the distribution P{f) for 3 and 
6 balls, which can be worked out analytically 19]. 

The two force networks shown in Fig. 1 are typical 
solutions / obtained by this scheme. We limit the dis- 
cussion to P{f) for interparticle forces, and address the 
boundary forces which show different distributions else- 
where 10, |23|. The interparticle P(/) for packings of 
increasing number of balls are presented in Fig. 2. Note 
that all P(/)'s display apeak for small /, which is typi- 
cal for jammed systems 'T"]. For large packings, this peak 
rapidly converges to its asymptotic limit. The tail of 
P{f) broadens with system size, but the present data is 
not conclusive about its asymptotic characteristics. 

Irregular packings We now apply our force ensem- 
ble approach to a more realistic system with a random 
packing geometry, and study a shear-stress induced un- 
jamming transition. To obtain a representative irregular 
contact geometry, we perform a standard Molecular Dy- 
namics simulation of a 50:50 binary mixture of 1024 par- 
ticles with size ratio 1.4 that have a purely repulsive 12-6 
Lennard- Jones interaction (a shifted potential with the 
attractive tail cut off)Lthe same system as the one stud- 
ied by O'Hern et al. [31 . We then quench such a finite 
temperature simulation onto a T = random packing 
with a steepest descent algorithm The static con- 

tact network that we obtain in this way then defines the 
matrix A in Eq. (Q. 

The P{f) obtained for this fixed packing is displayed 
in Fig. 121 even for a single contact geometry, we clearly 
reproduce a realistic P{f) which is very similar to both 
that of the triangular packings and to those obtained in 
experiments and simulations 0, Q . 



To investigate the role of the particle hardness, we have 
performed MD simulations of the same system with in- 
creasingly hard particles, obtained by varying the prefac- 
tor of the potential at constant pressure. For our original 
MD, which defined A, the particles are fairly soft with 
e ^ 0.1, and the corresponding P{f) is somewhat dif- 
ferent from the one in the force ensemble. When the 
hardness of the particles is increased, e diminishes and 
the corresponding P{f) indeed approaches the force en- 
semble P{f)- For the hardest particles (e ~ 0.02) these 
P(/)'s are virtually indistinguishable. We find that this 
holds for a variety of force laws [l9]. This confirms the 
validity of our approach for hard particles. 

Unjamming by shear It is also possible to study the 
effect of a shear stress on the force network ensemble by 
using the relation between the microscopic forces and the 
macroscopic stress field: 



(5) 



and extending the matrix of Eq. with the three linear 
constraints of Eq. . The average value of the force is 
set to unity by requiring a^x — cfyy = 1/2 HJ], and we 
vary r = cTxy/crxx- 

We find that P{f) evolves from a "jammed" distribu- 
tion with a peak, to an "unjammed" monotonous distri- 
bution as a function of shear stress (Fig. 4a). As a func- 
tion of the angle 0, (f ) varies in good approximation as 
l-|-2Tsin(20). This variation is consistent with Eq. jnj as 
well as with the alignment of the dominant contacts visi- 
ble in Fig. 4b - note the similarity to experimentally ob- 
tained sheared networks j^^. Since P{f) contains forces 
in all directions, the broadening of P{f) with shear stress 
follows immediately from this angular modulation of (f). 
However, this is only part of the story: we find that the 
shape of P{f) also varies with direction, from extremely 
jammed along the force lines, to almost purely exponen- 
tial along the weak principle direction (not shown) . 

When T approaches 1/2, (/) becomes zero along the 
weak principle direction, which implies that all forces 
along the weak direction approach zero. This can be 
interpreted as the breaking of contacts (see Fig. 4c). In 
addition, when t — > 1/2, the contact number drops to 
z = 4 (inset of Fig. 4a), and beyond this point, stable 
force networks are no longer possible. This simple mech- 
anism thus provides r = 1/2 as a definite upper bound 
for the critical yield stress P, 0, 0] ■ In terms of 
the "slip angle" this value corresponds to 30°. On the 
other hand, it has been speculated 0, that the qual- 
itative change of P{f) from peak to plateau could also 
be indicative of yielding; in our ensemble this occurs at 
T = 0.26 suggesting a slip angle of 15°. 

Outlook We have proposes a novel ensemble approach 
to athermal hard particle systems. The full set of me- 
chanical equilibrium constraints were incorporated, in 
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FIG. 4: Force networks under shear, (a) P{f) for increasing 
shear showing an "unjamming" transition at r « 0.26. The 
inset shows the contact number as function of shear, where 
contacts are considered broken when / < 10~*; for smaller 
values of the cutoff this curve remains essentially the same. 
(b) Examples of parts of the force networks under shear (c) 
Ratio of number of contacts with / > 10"'' and number of 
contacts as function of the contact angle show preferential 
breaking along the "weak" principal direction. 



contrast to more local approximations or force chain 
models 0, 0, I25I E^ . A number of crucial ques- 
tions can possibly be addressed within our framework. 
(1) Our approach is perfectly suited to include frictional 
forces, since these are difficult to express in a force law 
but simple to constrain by the Coulomb inequality. (2) 
The contact and force networks of sand piles exhibit dif- 
ferent anisotropics under different construction histories 
[isl E9I 1 . We suggest that contact network anisotropics 
may be sufficient to obtain the pressure dip under prop- 
erly created piles. (3) The problem defined by Eqs. Q 
could be generalized to arbitrary A and b, for which we 
can calculate P{f) and may ask under what conditions 
P{f) has an exponential tail, appear jammed etc. Pre- 
liminary work indicates that for realistic A but taking 
b nonzero with mean square average proportional to T 
captures the effect of a finite temperature of P{f) 
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